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Abstract 

The paper describes the aeroelastic analysis of a tiltrotor configuration. The 1/5 scale wind tunnel 
semispan model of the V-22 tiltrotor aircraft is considered. The analysis is performed by means of a 
multi-body code, based on an original formulation. The differential equilibrium problem is stated in 
terms of first order differential equations. The equilibrium equations of every rigid body are written, 
together with the definitions of the momenta. The bodies are connected by kinematic constraints, 
applied in form of Lagrangian multipliers. Deformable components are mainly modelled by means 
of beam elements, based on an original finite volume formulation. Multidisciplinar problems can 
be solved by adding user-defined differential equations. In the presented analysis the equations 
related to the control of the swash-plate of the model are considered. Advantages of a multi-body 
aeroelastic code over existing comprehensive rotorcraft codes include the exact modelling of the 
kinematics of the hub, the detailed modelling of the flexibility of critical hub components, and the 
possibility to simulate steady flight conditions as well as wind-up and maneuvers. The simulations 
described in the paper include: 1) the analysis of the aeroelastic stability, with particular regard to 
the proprotor /pylon instability that is peculiar to tiltrotors, 2) the determination of the dynamic 
behavior of the system and of the loads due to typical maneuvers, with particular regard to the 
conversion from helicopter to airplane mode, and 3) the stress evaluation in critical components, 
such as the pitch links and the conversion downstop spring. 
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Introduction 

The 1/5-scale aeroelastic model of the V-22 (Figure 1) was built and tested in the period from 1983 
to 1988, to support the preliminary design and the full scale development of the tiltrotor aircraft 
later known as the JVX. The wind tunnel tests began at the Transonic Dynamic Tunnel (TDT) 
of the NASA Langley Research Center (LaRC) on a semispan model, and were globally performed 
in three different facilities, including the Boeing Helicopter VSTOL tunnel for both the semispan 
and the full span model configurations. The semispan model is currently referred to as the Wing 
Rotor Aeroelastic Testing System (WRATS), and is located at the LaRC. 

The design of a complex system as a tilting rotor aircraft requires sophisticated analyses, with dif- 
ferent levels of approximation, as the design process matures from the preliminary, to the detailed 
analysis, up to the full scale development phase. Traditionally, rotors and rotor components are 
analysed with dedicated tools. Transfer matrix (a.k.a. Myklestad), FEA (NASTRAN, CAMRAD 
[11], [12], UMARC [14], DYMORE [1]) can be used for structural and modal analysis of the isolated 
blade; dedicated codes, based on simplified or imposed rotor kinematics, and relying on FEA for 
the modal characterisation of the flexibility of the rotor, can be used to find trimmed solutions of 
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the dynamics of the rotor. The family of codes that follow this philosophy is broad, and many of 
them have reached a high level of refinement and, expecially those developed in-house by helicopter 
companies, are successfully used in the design of commercial rotorcrafts (PASTA [15], DYN4-DYN5 
[23]). One further step towards generality is represented by comprehensive codes, that allow a more 
refined modelling of the kinematics as well as of the flexibility of the system, but with some sim- 
plifications, usually represented by assuming that the rotation speed is constant and the motion 
is periodic, leading to the direct solution of the trim problem (CAMRAD, UMARC), an iterative 
solution being usually required only when a wake model is used (CAMRAD, UMARC, CAMRAD 
II [13]). A unifying tool, able to cover such a wide spectrum of applications, is represented by a 
general purpose modelling code, based on a multi-body multidisciplinary approach, that gives the 
analyst modelling capabilities ranging over every design step without any undue simplification. 
Within this framework, one can progress from rigid body models, for preliminary performance 
definition, to fully detailed aero-servoelastic analyses, through intermediate steps, encompassing 
detailed mechanism definition, introduction of deformable elements, actuators, and control system 
components. This approach is likely to need more computer power than that required by spe- 
cialised and simplified approaches, but pays back in terms of efficiency since it allows the designer 
to avoid risky physical oversimplifications along with the greater modelling confidence allowed by 
using a single, general purpose, and well proven modelling tool. Moreover, with the computer 
power nowadays available, even the most complex models are likely to require a turnaround time 
that is compatible with an extensive set of parametric analyses and, within not so long a time, 
even with a complete system optimisation. Unfortunately, current commercial general purpose 
multi-body analysis codes, e.g. DADS [9], MECANO [3], ADAMS and others [25], still pose some 
limitations to the modelling of rotorcrafts, mainly due to insufficient aerodynamics, insufficient 
description of flexible bodies, and in some cases to limitations in the integration algorithms when 
applied to large finite rotations of the order of some revolutions [17]. 

The numerical simulations have been performed by means of an original multi-body formulation, 
which resulted in a code, MBDyn, that is still under development. It is intended for the simul- 
taneous solution of multi-disciplinar problems including non-linear dynamics, aero-servoelasticity, 
smart piezo-structural components, electric and hydraulic components. It allows the modelling of 
complex systems, a clear example of which is the proposed aeroservoelastic model of a tiltrotor 
aircraft. 

The multi-body formulation is presented first. The analytical treatment of the dynamics of rigid 
bodies, the handling of finite rotations, and the finite volume beam formulation are described, 
as well as some of the features of the multi-body code. Then the tiltrotor analytical model is 
described. The single components are detailed, and the steps of the build-up of the complete 
model are outlined. Results of intermediate analyses are presented, involving the kinematic and 
elastic characterisation of the hub and the modal analysis of the subparts. Finally some aeroelastic 
results, as well as the simulation of some typical maneuvres are presented. 


Multi-Body Formulation 

Dynamics. The multi-body problem is formulated by directly writing the core equations for each 
unknown (structural, electric, hydraulic, . . .). Constraints are imposed by means of appropriate 
equations, with “reacting forces” as unknowns, in a way that resembles the Lagrangian Multipliers 
method. It should be stressed that the reaction unkowns really are the reaction forces and not 
abstract multipliers. This solution has been preferred to a possible generalised variational approach, 
to simplify the writing of the equations and thus to reduce the computational burden. Moreover, 
a variational formulation for non-conservative loads and constraints, and for arbitrary elements 
in general, like those related to the previously mentioned non-structural elements, would have 
been too cumbersome. For computational efficiency and simplicity of formulation, the dynamics 
problem has been stated as a first order differential system of equations; the equilibrium equations 
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of a body are: 

j q — F ( x , x, R, to, . . . j t) =0 
\ 7 — (to x S) x x — M (x, x, R, u >, , . . , t) =0 

where q, 7 are the momenta, x is the position of the node, R is the rotation matrix that describes 
the rigid rotation from the local to the global frame, to is the angular velocity of the node, that is 
related to the rotation matrix by the relation ux = RR T ' 2 . The inertial forces and moments balance 
in a d’Alembert sense the forces F and moments M, which may depend on the configuration (i.e. 
on the kinematic unknowns, to give elastic forces) as well as on other parameters (like the time 
t). The term (cu x S) x x in the moment equation is required since moments are referred to the 
position of the node, that is variable. The definitions of the momenta must be added, as follows: 

j mx + S x lu — q = 0 
\ —S x x + Juj — 7 = 0 

The mass of the body is m; the inertial properties S, J represent the first and second order inertia 
moments of the body and are referred to the global frame. Their transformation from the local 
to the global frame is given by S = RS, J = RJR T . The physical momenta have been chosen 
as unknowns to let different elements with mass contribute to the inertia of a node. Kinematic 
constraints are added as constraint equations $ {x,x: t R,uj, . . . ,t) = 0. The unknown constraint 
reactions Vp, Vm then contribute to the equilibrium equations as Lagrange multipliers. The system 
is summarised as follows, in a form known as Lagrangian of the first kind: 

7nx + S x to — q = 0 
—S x x + Jlo — 7 = 0 
< q + Vp — F (x, x, R, u ), . . . , t) =0 

7 — (uj x S) x x + Vm — M (x, x,R,lo,. . . ,t) = 0 
•l> (jyj\ />’. _c. . . . . /) (I 

The handling of the rotation matrix R is detailed in the following section and in Appendix A; it will 
be expressed in terms of (arbitrary) rotation parameters g. Equations $ represent holonomic and 
non-holonomic constraints. In general these equations are algebraic or mixed differential-algebraic 
for holonomic and non-holonomic constraints respectively, and thus the system of equations is 
Differential Algebraic (DAE) of index three [2], It is important to remark that the system is solved 
“as is”, in the unknowns x, g, q, 7 , V>, Vm without any further substitution. Initial satisfaction of 
the constraint equations is ensured by properly assembling the joint elements in a system in which 
virtual springs link the nodes to their assigned positions, and then iterating until a compatible 
solution is obtained. This allows a great flexibility, since the initial (given) configuration may be 
non-compatible. Eventually the system is trimmed in a configuration that satisfies equilibrium by 
solving it at the initial time, with a modified update procedure that preserves the configuration and 
modifies only the momenta and the reaction unknowns. In this way initial equilibrium can be com- 
puted without differentiating the constraint equations up to the second order. Finally the solution 
is improved by iterating the so called “fictitious steps” , that are normal solution steps performed 
at a very small time step and with high numerical damping, followed by a single-step restart of the 
computation. The fictitious steps implicitly perform a sort of numerical differentiation of the con- 
straint equations without resorting to any specific procedure. General elements, called Genel, are 
provided to account for user-defined generic components. They introduce user-defined unknowns 
a, called abstract degrees of freedom, that can interact with standard pre-existing unknowns. As 
a result, forces and moments may depend on abstract unknowns too, and Genel equations are 
added to the system in the form: 

j F, M = F, M (;c, x,R,ui, a, a , . . . , f) 

| G ( x , x , R, lu, a, a , . . . , t) =0 

They are used in the preliminary development phase, to allow a high degree of flexibility in writing 
the formulation and the code; when a reasonable standardisation is reached, they are merged 

2 Operator (■) x allows to write the vector product in matrix notation, i.e. being a a vector, ax is the matrix 
that multiplies vector b to give a x b 
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in a separate family of elements. The system is being extended to the integrated simultaneous 
solution of multi-field problems; the active control of smart structures has been already addressed 
in [5] by adding electric variables and sensors, actuators and signal processing elements. Further 
development will involve the reinclusion of modal coordinates, and the modelling of hydraulic 
circuits and elements, to simulate the control system of a complete rotorcraft. 

Kinematics. The formulation of the multi-body problem is based on the kinematic description 
of large displacements and rotations. For an accurate but efficient, i.e. fast and cheap, solution, 
an appropriate formulation for the rotations is fundamental. In this work, the total position of 
each body is used as unknown. The choice of the rotation unknowns pose a delicate problem. 
The orientation of a local frame is described by an orthogonal matrix R that maps vectors from 
the local to the global frame. The parametrisation of large rigid rotations requires at least three 
unknowns, but four parameters are needed to avoid singularities that arise when the orientation of 
the rotation is undefined. This problem has been prevented by considering incremental rotations as 
unknowns. So the current orientation of a reference frame R is accounted for by a constant rotation 
matrix R r , that multiplies an unknown incremental rotation, held by matrix R&, that is assumed 
to be small enough to avoid any singularity. Matrix R r is updated after every step. A widely used 
three parameter parametrisation is represented by the rotational vector ip, that describes a finite 
rotation of amplitude given by its modulus, about an axis represented by the vector itself. A very 
efficient parametrisation of the finite rotations has been found in the Gibbs-Rodriguez parameters 
<7 = 2 tan -1 (;p/2), modified with respect to the conventional notation by means of the factor 2. 
This makes the linearised expressions of the rotational entities concide with those of the rotational 
vector Lp. The rotation matrix R is: 


R = r + , , 4 T (fl X X <?x ) 

4 + g 1 <? V 2 J 

Gibbs-Rodriguez parameters do not use trigonometric functions, that are computationally expen- 
sive. Some useful properties of the rotation matrix, and its differentiation rules, are reported in 
Appendix A. 

Integrator. Time integration is performed by means of an implicit, A/T-stable, second order 
accurate, predictor-corrector integrator. The basic formulas are: 

12 

Vk = — ^ (Vk-1 ~ Vk-i) + 8 yk-i + 2 

Vk = (i - +Pyk-2 

+ h ^<5 + jjk + h ^-/ 3 + - - 2<5^ jf*_i + h ^-/3 + <5^ yu-2 

for the predictor, consisting in the cubic extrapolation of the derivatives based on the states 
and their derivatives at the two preceding steps, and the prediction of the state based on an 
extrapolation, the coefficients of which ensure second order accuracy, with user-defined control of 
the numerical damping; h is the time step. The formulas have been generalised to a variable step 
predictor; in this work, for sake of simplicity, only constant time step formulas are described. The 
coefficients f3 and S can be expressed in terms of the desired asymptotic spectral radius , under 
the assumption of real and coincident asymptotic roots (best fit criterium): 

n = 4 f4 ~ (1 "Poo) 2 c = (1 ~Poc) 2 

4 — (1 — Poo) 2 2(4— (1— Poc ) 2 ) 

For poo = 1 the method is very similar to the Crank-Nicholson rule (no dissipation), though 
using two steps, while for p 0 0 = 0 the method coincides with the well known Backward Difference 
Formulas (BDF) [2], The correction is performed by means of a complete/modified Newton- 
Raphson iteration. 

Rotations Updating. The correction is made referring to the predicted reference (this technique 
has been called updated-updated). The predicted frame is used as reference, so it is constant; only 
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the rotation related to the correction is unknown. As a consequence, the unknown rotation becomes 
really small, provided the prediction is accurate. The angular velocity can be expressed as: 

L0 = G'a<7A T R&iOr 

where subscripts (-) A mean that the rotation parameters and their derivatives are referred only 
to the correction rotation. As a consequence, when the prediction is accurate, the terms involved 
in the linearisations can be approximated as: Ga — I, Ra — I, AGa — 0. The differentiations 
become: 

A R = Ag A x R r , A to = Ag A + Ag A x to r 

this greatly simplifies the writing of the Jacobian matrix, with consequent savings in computa- 
tional time. Accuracy is preserved by consistently calculating the residual. This approach can 
be regarded as an “intrinsically modified” Newton-Raphson solution, since the Jacobian matrix is 
approximated during the construction, and thus acts as a sort of very well suited preconditioner. 
Constraints. Many constraints can be formulated. The most commonly used have already been 
implemented, to allow the modelling of common mechanisms. Table 1 lists the available elements. 
Most of them are based on the orthogonality condition of two vectors, namely e' a e J b = 0, where 
the subscripts refer to the bodies the vectors belong to, and the superscripts refer to the directions 
in the local reference frame. In Appendix B, the linearisation of basic constraints is described. 
Beams. Beams represent a fundamental elastic element in the proposed multi-body implemen- 
tation. The kinematic description of the generalised deformations, i.e. strains and curvatures, is 
based on an intrinsic formulation of the beam. The strains are defined as the difference between 
the current and the initial derivatives of the reference line p (£) that describes the position of the 
beam. The strains, referred to the material frame, are: 

e = R 0 R T p' - p'o 

where the position p refers to the current frame, while po refers to the initial configuration of the 
beam, the prime (•/ performs a spatial derivative, and Rq is the rotation matrix at point £ in the 
initial configuration. The linearisation of the strains results in: 

Ae = RqRJ (A p' + p' x Ag A ) 

The geometric curvature p is defined as the spatial derivative of the reference frame of the beam 
section: 

px = R t R' 

The difference between the current and the initial, or imposed, curvature po, represents the elastic 
curvature ft , expressed in the material frame: 

ft x = R T R' - RqR'q = p x -pox 

When incremental rotations are considered, the elastic curvature becomes: 

ft = {R T &g T R^r 

After linearisation: 

Ak^RoR^Ag' A 

The internal forces are defined by means of a constitutive law in terms of generalised strains, i.e. 
•d = •& {ip) where d are the internal forces and moments and ip are the strains and the curvatures. 
An original finite volume approach is used to formulate the beam elements. A linear elastic 
constitutive law is usually considered, but, without any loss in generality, an arbitrary constitutive 
law can be considered as well, since the finite volumes formulation deals with collocated internal 
forces only, irrespective of their nature. Finite volumes applied to the equilibrium of a beam can 
be interpreted, in physical terms, as the direct balance of the forces that act on a finite piece of 
beam. The equilibrium equation of a finite piece of beam is: 

( I - U b ) i) b - (I - U a ) i) a = T b a 
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where a and b label the ends of the piece of beam, T h a are the resulting dead forces and moments ap- 
plied in the interval [a, b\, and matrix U represents the arm of the internal forces in the equilibrium 
equation of the moments: 


u(0 


0 0 
- ( V (0 -x)x o 


being x the pole the moments are referred to. A three node parabolic C° element has been 
implemented; it gives the exact solution for end-applied loads [18], [6]. The inertia of the beam 
is accounted for in a lumped scheme. Consistent inertia forces have been formulated, and their 
implementation in the linearized case showed higher accuracy in modelling the dynamics of a 
coarse beam model, but at a higher computational cost [6]. The finite volume beam can be 
regarded as a constraint that relates the reaction forces to the deformation of the link, and thus 
to the configuration of the system. Provided the relation between reactions and configuration 
is invertible (i.e. the Hessian of the strain energy is positive definite), the constraint equation 
can be implicitly solved, thus allowing the direct writing of the contribution of the beam to the 
equilibrium equations in terms of position and rotation unknowns. The loss in symmetry of the 
linearised matrices of the finite volumes beam is not a drawback in the context of a multi-body 
formulation, since the problems at hand already are non-symmetric. Finite volume beams are easy 
to implement in a multi-body formulation since only collocated evaluation of contributions to the 
equilibrium equations is required. Moreover, they straightforwardly resemble the natural partition 
in distinct bodies that is peculiar to the multi-body formulation. The finite volume description of 
the deformation of slender bodies is consistent with the mathematical, intrinsically discrete, multi- 
body model, and thus allows an easy but thorough modelling, especially when used in conjunction 
with specialised beam section analysis and characterisation formulations, like those described in 

[7] and [4], 

Aerodynamic Forces. The aerodynamic forces are based on the strip theory, using elements 
that refer to rigid or beam shaped aerodynamic surfaces. The aerodynamic coefficients are based 
on the interpolation of experimental data spanning 360 degrees of angle of attack. Corrections 
are made to determine the drag due to spanwise flow, as well as the effects of dynamic stall 

[8] . Quasi-steadv aerodynamics coefficients can be considered as an alternative approach. Rotor 
elements are defined, to account for the effect of rotor induced velocity with an increasing degree 
of refinement, from uniform up to Glauert and Mangier inflow distribution in hover and in forward 
flight, respectively. Dynamic inflow modelling can also be used [22], The previous version of the 
code was interfaced to an aerodynamic code that models the wake of a rotor. An upgrade for the 
current version has been planned. 


Tiltrotor Models 

A tiltrotor aircraft is a complex system that has the behaviour of both a conventional airplane and 
of a rotorcraft, with peculiar maneuvres, e.g. the conversion. Its aerodynamics are characterised 
by the high influence of the rotor on the airstream that affects the wing, in both helicopter and 
airplane configuration. In fact the WRATS project has been mainly focused on the reduction of the 
vibration level induced by these interactions in the airplane mode by means of an actively controlled 
swashplate with the Higher Harmonic Control (HHC) technique, coupled to an active flap [20]. The 
blades of the rotor represent a compromise between helicopter and propeller blades. Since they are 
optimised for the high axial airstream speed typical of the airplane mode, they are highly twisted, 
thus showing high elastic couplings between twist and in- and out-of-plane bending [19]. The 
conversion maneuvre, when performed at typical rates for aircraft control, introduces gyroscopic 
effects that are unusual in conventional helicopters. The flexibility of the wing can magnify the 
effects of ground and air resonances, the latter being typical of automatically controlled rotorcrafts. 
Many of these problems are well known, but have never been faced to this extent, while others 
are completely new. For this reason, a particularly intensive experimental campaign preceeded 
and accompanied the whole development of the JVX [16], supported by numerical analyses of 
the related subproblems. Eigenvalue analyses of the rotor dynamics, the determination of flutter 
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margins of the rotor, of the wing and of the ensemble, by means of analytical models based 
on comprehensive rotorcraft analysis codes, and dynamic simulations of the rotor mechanisms by 
means of early multi-body codes have been performed at each step of the wind tunnel investigations 
[23], [24], Most of these analyses were performed by means of comprehensive rotorcraft codes, with 
properly simplified models. In this paper a different approach is proposed, consisting in the use 
of a single, state-of-the-art, multipurpose multi-body code that is able to perform most of the 
required simulations starting from a single bulk model that can be specialised for each analysis. 
Each subpart of the tiltrotor is modelled and analysed in its basic kinematic and dynamic features, 
then the parts are assembled together and the system is analysed as a whole. By using the same 
code and the same modelling for the single parts and the assembly, and by using a rather general 
approach in the kinematic and mechanical description of the single parts, any undue approximation 
is avoided. The tiltrotor has been split in the following subsystems: 

1. The blade, made up of a flexbeam, a pitch hinge and a pitch link. Either rigid or flexible 
blades have been considered. 

2. The gimbal, a constant velocity joint made up of its components in order to give an accurate 
kinematic description of the joint. 

3. The swash plate, made of the two plates, the two scissors that constrain the axial rotation 
of the plates with respect to the pylon and to the hub, and the three non-rotating links that 
control the collective and the cyclic position of the plates. 

4. The half-wing model, made of the deformable wing, the pylon, modelled as a rigid body, the 
conversion hinge, the downstop spring and the mast. 

The complete model is sketched in Figure 2. 

Blade. The single blade model has been used to analyse the dynamic properties of the iso- 
lated blade, such as frequencies and aerodynamic properties. Three different models have been 
considered, with increasing discretisation refinement. All the models share the description of the 
flexbeam, that uses a three-node beam element, and of the controls. The blade is joined to the 
flexbeam by a spherical hinge at the outer end, and by a spanwise oriented in-line joint, 2.2 in. 
outwards of the rotor axis. The bending of the flexbeam accounts for flexible flap and lag motion, 
while the pitch rotation is allowed by the bearings. A distance joint between the rotating swash- 
plate and an offset point aft of every blade cuff models the control link. It can be both rigid or 
flexible, to account for the flexibility of the control system. A rigid blade has been considered first, 
which proved to be poor because the yoke is very stiff, so even the very first in-plane bending mode 
implies appreciable deformation of the blade itself. On the other hand, the rigid blade represents 
an acceptable tradeoff when the performances of the rotor are addressed. A flexible blade has been 
subsequently modelled, with two and four beam elements respectively. The first frequencies of the 
cantilevered blade obtained by MBDvn are reported in Table 2, compared to Ground Vibration 
Tests (GVT), UMARC and NASTRAN Finite Element Analysis codes results. 

Gimbal. The gimbal model has been used to determine the kinematic and gyroscopic properties 
of the rotor. It consists in a constant velocity joint, made by two universal joints, linked to the mast 
and to the hub respectively at one arm of each cross. The other arm of the crosses is connected to 
a linkage, that transmits the torque between the mast and the hub and keeps constant the distance 
between the two U-joints. The hub is also linked, by means of an in-line joint, to a spherical joint 
on the mast that allows the gimbal motion. The gimbal allows the rigid flapping of the whole 
rotor and, since the direction of the angular velocity tilts together with the hub, no Coriolis forces 
due to this motion result in the blades when the flapping is steady. At the same time, the 1 per 
rev. flapping motion has no stiffness due to centrifugal effects, but only that provided by a set of 
springs. 

Swash Plate. The swash plate model has been used to analyse the kinematics of the control 
system and, together with the gimbal and the rigid blade, to evaluate the pitch-flap-lag couplings 
for the whole collective pitch range. It has been used also to apply the desired controls to the 
rotor during dynamic simulations. It consists in the two plates, modelled as rigid bodies joined by 


7 



a plane hinge. The fixed plate is linked to the pylon by means of an in-line joint that forces it to 
translate along the mast. A fixed and a rotating scissor constrain the rotation of the two plates 
about the rotor axis, with respect to the pylon and the mast. Three variable distance joints are 
used to control its translation (collective pitch) and attitude (cyclic pitch). The elongation of the 
fixed control links is imposed by means of a dedicated general purpose element that splits the three 
fundamental control inputs, namely collective, and fore/aft and lateral cyclic pitch angles, into the 
elongations of the links. Figures 3-5 respectively represent the kinematic pitch-flap coupling due 
to the gimbal and the flexible flapping, and the control stiffness as function of the collective pitch, 
compared to data obtained from Bell Helicopter and model calibrations. 

Wing-Pylon. The half wing model has been used for aeroelastic clearance of the isolated 
wing. Both dynamic and aeroelastic properties of the wing in forward flight have been analysed. 
It consists in two beam elements for the wing, and in the pylon, modelled with a rigid body. The 
pylon is connected to the wing by means of a flexible spindle, modelled with a beam element, and 
a downstop spring. The spindle models the conversion bearing. Its bending allows the pylon to 
rotate with respect to the wing about the roll and yaw axes. The conversion actuator constrains 
this rotation, and controls the conversion angle. In the wind tunnel model, springs with different 
properties are used to simulate the behavior of the conversion actuator in helicopter and airplane 
configuration, both on- and off-downstop. The main frequencies of the wing-pylon are reported in 
Table 7, compared to both GVT and NASTRAN results. 

Some preliminary considerations on the flapping motion: the rotation center for the flapping 
due to the gimbal, namely 1 per rev. flapping, is located on the rotation axis, 2 in. above the 
rotor plane, while the one for the flapping due to the flexbeam deformation, namely the cone 
and >1 per rev. flapping, lies about 1.5 in. outboard along the blade axis. The pitch control is 
linked to the blade 75° aft of the blade itself and thus introduces a pitch-flap coupling S 3 = —15° 
that is negative (flap up causes pitch up) for the 1 per rev. motion, and slightly positive for the 
flexible flapping. It is known, [10], that, for a stiff-in-plane rotor, a positive S 3 can give raise to a 
pitch-flap aeroelastic instability when the first out-of-plane and in-plane frequencies nearly meet. 
The occurrence of this instability in the simulations required a deeper analysis of the flexibility 
of the yoke. In detail, the flexibility of the inner part of the yoke, from the hub to the inner 
pitch bearing, proved to be fundamental in describing the correct coupling between the blade pitch 
and the flexible, symmetric flapping motion. After this part was properly modelled, the slight, 
symmetric instability in the analytical model disappeared. Figure 4 shows the change in pitch-flap 
coupling for the cone flapping motion both for rigid and flexible root of the flexbeam. 

Rotor Models. The model, consisting in the rotor with rigid/deformable blades, the gimbal 
and the swashplate, has been used to investigate the stability of the rotor, with particular regard 
to the pitch-flap-lag coupling, and to evaluate the aerodynamic response to the controls. The 
rigid blade model matched the out-of-plane frequencies, but gave poor results for the in-plane 
frequencies, so it was of little use in other than performance analyses. The flexible blade models 
agreed very well with available data for the low frequencies of the rotor, both in the rotating and 
non-rotating cases in vacuo. Both the single blade and the complete, three blade rotor models have 
been analysed. For this purpose, UMARC has been modified to allow the modelling of multiple 
blade rotors in the finite elements analysis module, with multiple load paths to account for the 
control links. Tables 3 and 4 show that the presence of the gimballed hub modifies the natural 
frequencies of the system by breaking the symmetry. In fact, non-symmetric modes are found, 
as shown both by the analysis and the experiment. Results from Bell Helicopters were available 
for the locked gimbal case, since they were obtained for a single blade model. They refer to an 
old configuration of the hub, with calibrated springs at the blade root to simulate the stiffness of 
the controls. These results are not completely representative of the current configuration. The 
GVT results with locked gimbal are also not completely significant, since the gimbal couldn’t be 
perfectly locked. As an example, in Tables 5 and 6, the rotating frequencies of the complete, three 
blade rotor are reported at two typical rotating speeds and collective pitches, respectively referring 
to hover and cruise flight conditions. 

Wing Model. The wing model shows good correlation for the lowest modes, as reported in 
Table 7. The results are compared to experimental measurements and to numerical results based 



on NASTRAN code [21]. The beam and torsion inodes are strongly coupled and are influenced by 
the properties of the downstop, the spring that is used to lock the conversion actuator in airplane 
mode. At present there is no conversion actuator on the wind tunnel model, so it is simulated by 
a set of springs that model its stiffness in different configurations. 

Wing-Rotor Models. The previously mentioned models have been merged by mounting the 
rotor models on the flexible wing. In the rigid blade version, the complete model has been used to 
evaluate the performances of the aircraft during maneuvres, an example of which is the conversion. 
The deformable blade model has been used to test the stability of the elastically mounted rotor and 
to assess the feasibility of the multi-body model for the simulation of the whole, detailed deformable 
system. Moreover, the effects of the flexibility of the blades on the dynamics of the system, in terms 
of transmission of the higher harmonics of the rotor system to the body of the aircraft, have been 
investigated. The first wing modes are not directly affected by the modelling of the flexibility of the 
rotor. The torsion mode of the wing is very close, and at some airstream speeds coincident, to the 
rotor speed; this gives raise to resonance that can be seen in the frequency analyses of the internal 
forces of the wing. Four wing modes are mainly considered: the beam, chord and torsion modes of 
the wing, and the so called “pylon yaw” mode, a low frequency yaw oscillation of the pylon due to 
the flexibility of the conversion actuator. When considered in the fixed frame, the retreating rotor 
modes interact with the wing modes. This can be clearly appreciated from a frequency analysis 
of the wing response when the rotor modes are excited. Most of these modes cannot be easily 
identified when the aerodynamics are modelled, since they are highly damped. For this reason, a 
comprehensive analysis of the structural properties of the model has been performed by simulating 
in vacuo operations, while the aeroelastic properties have been estimated in different ways. The 
damping of the wing modes in forward flight has been estimated by system identification of the 
(damped) response to a given input, as is usually done during actual wind tunnel tests, while 
the aeroelastic pitch-flap coupling has been estimated by measuring the phase shift between an 
harmonic control input and the flapping response. 

Figures 6, 7 refer to a collective pitch maneuvre. They show the internal moments at the wing 
root and the geometric pitch of blade #1 as the collective control is raised from 0 to 10 degrees in 
one second. The simulation is performed in helicopter mode; the nominal hover rotation speed of 
888 rpm is reached in one second to obtain a trimmed condition (not shown). There is no airstream 
speed. The difference between the given control and the actual pitch of the blade is due to the 
deformation of the flexbeam and of the flexible link. 

Figures 8, 9 refer to a 5 degrees fore/aft cyclic pitch maneuvre. As the rotor tilts forwards, 
the high frequency in-plane modes of the wing are excited, as shown by the plot of the internal 
moments. The oscillations in the pitch link are 1/rev., partially due to 1/rev. flexbeam flapping 
that is superimposed to the gimbal flapping (which implies no appreciable pitch link loads), that 
is needed to counteract the gimbal springs. 

Figure 10 refers to the conversion maneuvre performed by a deformable blade model. It shows 
the internal forces at the wing root. The conversion is performed at a 10 deg/s constant angular 
speed. Oscillations of the internal forces due to the untrimmed initial conditions are appreciably 
damped as the maneuvre proceeds to the end, at 9 s. The following abrupt raise of the internal 
forces is due to the transient caused by the end of the maneuvre. 

The flexible blade model has been used to simulate the response to a cosinusoidal vertical gust in 
airplane mode, of 10 ft/s amplitude. Both the stability and the sensitivity of the tiltrotor have been 
addressed. Figures 11, 12 show the wing out-of-plane internal moment due to the gust at different 
airstream speeds, for both off- and on-downstop configurations. In figure 11 the off-downstop 
configuration is clearly approaching the stability boundary at 140 Kts. When the rotating speed is 
increased, the stability boundary moves towards lower speeds, as shown by previous analyses and 
experiments [23]. 

The deformable blade model has been used in the longest and most sophisticated performed 
simulations, i.e. the complete conversion. The model is made of 45 nodes, 39 rigid bodies, 31 
joints of different kind, 18 beam elements, 14 aerodynamic elements for a total of 571 degrees of 
freedom. Initially, the rigid blade model required At = .5 x 10 -3 s, while the deformable blade 
model required At = .25 x 10 -3 s to start correctly. When a variable time step was used, during 
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the rigid blade model simulations it grew very quickly to 3.0 -t- 3.5 x 10 -3 s, while, during the 
deformable blade model ones, it reached about 1.0 -t- 1.2 x 10 -3 s. The conversion simulation 
required about 4.5 hours on a Pentium PRO 200 PC for a total of 40000 fixed size time steps (10 s 
at At = .25 x 10 -3 s). When performed with variable step size, it required about one hour. After 
the model has been refined, and a soft start has been used, the flexible blade model is able to start 
with At = 10 -3 s, requiring 1.1 hours on a Pentium II 333 PC. 


Concluding Remarks 

An original formulation for the solution of the multi-body multidisciplinary modelling of a tiltrotor 
has been presented. The formulation proved to be efficient without excessive simplifications. Effi- 
ciency has been enhanced while maintaining the physical meaning of both the equations and the 
unknowns. On the computational side, the redundancy of the formulation is exploited by a sparse 
solver, while the updated-updated rotation approach both reduced the computational burden and 
eased the writing of the jacobian matrix required by the nonlinear solver. Correctness is achieved 
by avoiding any undue approximation in writing the kinematic relations and by properly calculat- 
ing the residuals even in presence of computationally convenient approximations of the Jacobian 
matrix. The use of an unconditionally stable method allows a time integration whose step is dic- 
tated essentially by the desired accuracy. A model of the tiltrotor used in WRATS investigation has 
been analysed, consisting in rotor models of increasing refinement, with rigid and flexible blades, 
the gimballed constant velocity joint, the swashplate and the control links, the pylon, the conver- 
sion hinge and the flexible wing. Aerodynamic loads have been considered, to simulate different 
test conditions, from aeroelastic stability investigations to the simulation of complex manoeuvers, 
including conversion and blade pitch control. The adopted approach represents an interesting way 
to perform analyses of complex systems, of which the tiltrotor model is a clear example, and can 
lead to a valid and thorough design tool. It is worth underlining that the results here presented are 
meant to assess the validity of the formulation and of the code, rather than a complete analysis of 
the WRATS tiltrotor model. The simulations have been performed on a Pentium Pro 200 and on 
a Pentium II 333 PCs using Linux OS. In view of the proof-of-concept code perspective, the basic 
problem solving engine is retained to be significant a competitor with respect to corresponding 
state-of-the-art commercial codes, but it lacks in the pre- and post-processing features of such 
softwares. Since, even if substantially sound and robust, Linux compilers are not as effective as 
proprietary compilers used on the most powerful available RISC workstations, it is expected that 
significantly better performances can be obtained in a real industrial design environment even 
with the code “as is”. Nonetheless enhancements are planned for the adaptive time step control 
in the integration, for more effective iterative solvers and for a rough low cost coarse scale parallel 
implementation on PCs using message passing paradigms. 
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Appendix A: Finite Rotations 


The differentiation of a vector v = Rv that is constant in a local frame is: 

v' = R'v = R'R t v 


where the orthogonality of matrix R has been exploited, that is, the scalar product of two vectors 
does not depend on the reference frame: a T b = a T b , from which follows that R T R = I. When 
total Gibbs-Rodriguez rotation parameters g are used, the following expression holds: 

R'R t = ( Gg ') x 


with matrix G given by: 

G = iv7s( /+ i" x ) 

Since matrix G depends on the rotation parameters only, the derivative of matrix R linearly de- 
pends on the derivatives of the rotation parameters. When differentiating with respect to time, or 
linearizing, g or A g must be respectively used instead of g 1 . Then the differentiation of matrix R 
assumes the meaning of angular velocity u = Gg and of rotation perturbation 6 a = (GAg) respec- 
tively. It is apparent that if parameters g are assumed to be small, matrix G can be approximated 
with an identity matrix /, in the spirit of the updated-updated rotation approach. 


Appendix B: Simple Constraints 

As an example, two basic constraints are outlined. Most of the constraints can be obtained as a 
combination of the coincidence and orthogonality conditions here presented. 

1 . Coincidence. Let Xi and p i = 1 , 2 , represent the position of two independent nodes 1 , 2 
and the offset from the nodes to the position of the joint, both in the global frame. For sake 
of simplicity, the offsets pi are assumed to be constant in the local frame. The constraint 
equation is: 

(x-2 + p2 ) - (X! + pi) = 0 

This constraint generates a reacting force / at the coincidence point, which on turn results 
in a force at the two nodes, and a couple due to the offsets: 

F = ±f M = ± Pi x / 

The differentiation of the constraint and of the reactions gives: 

(Ax- 2 - p2 x 6A2) - (Axi - pi x 6A1) = - (X2 + P2) + {xi + pi) 

±A / = =F / ± (Pi x A / - / x pi x 6 a i) = FPi x / 

where 6a refers to the perturbation of rotation, namely GAg. 

2 . Orthogonality. Let e,, i = 1 , 2 , represent the unit vectors of some coordinate direction 
referred to two independent nodes 1,2 and expressed in the global frame. The orthogonality 
constraint equation is: 

62 ei = 0 

This constraint generates a reacting couple to that is scalar, and acts in direction eo xei: 

M = ± (eo x ei) to 

The differentiation gives: 

(e 2 x ei) T #A2 - (e 2 x e\) T 6 ai = -e-fei 
A to (ei x e 2 x 6a2 — e 2 x e\ x #ai) ± (e 2 x e\) Am = T (e 2 x ei) to 

By adding 1 , 2 and 3 constraints of this kind, universal, plane and prismatic hinges can be 
respectively obtained. 
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Table 1: List of available elements 


Aerodynamic 

Rigid 

Beam 

Induced Velocity - Rotor 

Dynamic Constraints 

Beam 

Rod 

Spring 

Hinge 

External 

Force 

Couple 

General 

Generic Sensors Dynamics 
Swash Plate Controls 
Discrete Time Linear Control 
User Defined Element 

Inertia 

Lumped Mass 

Kinematic Constraints 

Clamp 

Distance 

Spherical Hinge 

Universal Hinge 

Plane Hinge 

In-Plane 

Prismatic 

Imposed Velocity 

Imposed Angular Velocity 
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Table 2: Cantilevered blade frequencies , Hz 


Mode 

Exp 

UMARC 

NASTRAN 

MBDvn 

4 elem. 2 elem. 

1 Beam 

12.29 

12.3 

11.5 

11.3 

11.7 

1 Chord 

34.11 

34.1 

33.4 

33.1 

32.7 

2 Beam 

52.44 

53.0 

56.7 

55.8 

55.0 

1 Tors. 

113.35 

111.4 

127.0 

119.0 

122.0 
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Table 3: Single blade with flexbeam (locked gimbal), non-rotating, Hz 


Mode 

BELL 

GVT 

UMARC 

MBDyn 


10 deg. 

50 deg. 

10 deg. 

10 deg. 

10 deg. 

50 deg. 

Cone 

6.6 

7.2 

6.6 

6.3 

6.8 

7.8 

2 F 

26.6 

35.1 

25.2 

30.8 

28.5 

39.0 

3 F 

68.8 

77.3 

69.3 

77.9 

73.5 

82.0 

1 L 

19.3 

12.5 

20.6 

19.3 

19.5 

12.6 

1 T 

114.5 

109.9 

112.6 

110.0 

109.0 

107.0 
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Table 4: Full rotor (free gimbal), non-rotating , Hz 


Mode 

GVT 

MBDyn 


10 deg. 

10 deg. 

50 deg. 

Gimbal 

2.0 

1.8 

1.5 

Cone 

6.8 

7.0 

7.8 

2 Flap 

25.0 

26.5 

36.5 

2 Flap asym. 

64.2 

57.1 

55.0 

3 Flap 

76.2 

78.0 

82.5 

1 Lag 

19.7 

19.0 

12.7 

2 Lag 

91.3 

98.0 

92.0 

1 Torsion 

112.1 

109.0 

107.5 
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Table 5: Rotating frequencies, 888 rpm, 6 75 % = —3 deg, Hz 


Mode 

Myklestad 

UMARC 

MBDyn 

Gimbal 

- 

14.8 

14.8 

Cone 

17.2 

17.3 

17.5 

1 Lag 

22.4 

20.8 

24.0 

Coll Lag 

42. 

44.0 

36.0 

2 Flap 

37.33 

49.6 

41.0 

2 Flap asym. 

- 

70.2 

65.0 

3 Flap 

75.33 

90.3 

73.0 

Flap/Torsion 

89.33 

92.7 

mm 

Lag/Torsion 

- 

113.4 

rnBBm 

Torsion 

- 

116.0 

HB9 
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Table 6: Rotating frequencies, 742 rpm, #75% = 55 deg, Hz 


Mode 

Myklestad 

UMARC 

MBDyn 

Gimbal 

- 

12.4 

12.6 

Cone 

14.7 

14.9 

15.1 

1 Lag 

15.3 

15.8 

16.5 

2 Flap asym. 

- 

42.3 

44.2 

Coll Lag 

32.7 

45.9 

46.9 

2 Flap 

45.3 

45.6 

49.1 

3 Flap asym. 

- 

46.9 

60.3 

3 Flap 

66.0 

60.1 

65.2 

Flap/Torsion 

89.3 

90.6 

97.8 

Lag /Torsion 

90.0 

90.8 

89.7 

3 Lag 

- 

92.0 

92.9 

Torsion 

- 

116.0 

108.5 
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Table 7: Wing frequencies, Hz 


Mode 

downstop 

GVT 

NASTRAN 

MBDvn 

Beam 

-'ESI 

6.00 

6.16 

5.9 


.El 

5.51 

5.45 

5.4 

Chord 

1*9 

8.45 

9.33 

9.1 


off 

8.45 

8.74 

0° 

bo 

Torsion 

on 

12.5 

12.6 

12.5 


off 

10.6 

10.6 

11.0 

Pylon Yaw 

on 

16.5 

18.9 

17.2 


off 

16.7 

16.7 

16.6 


19 










Figure 1: WR ATS Model at Langley’s TDT 
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Figure 2: Analytical Model 
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Collective Pitch at 75% of the Blade, deg 


Figure 4 : Pitch-cone coupling as function of the collective pitch #75% 
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Control Stiffness, in-lb/rad 



Figure 5 : Control stiffness as function of the collective pitch #75% 
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Figure 11: Out-of-plane bending moment of the wing due 
to a vertical gust of 10 ft/s at airstream speeds 100-lfO 
Kts , flexible blade model 
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Figure 12: Out-of-plane bending moment of the wing due 
to a vertical gust of 10 ft/s at airstream speeds 160-200 
Kts , flexible blade model 





